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Introduction 


Users  of  least-squares  collocation  ask  for  a theory 
that  gives  an  answer  to  practically  meaningful  questions:  What 
is  the  accuracy  of  our  results?  Can  we  apply  statistical  testing 
techniques?  How  can  we  compute  statistical  distributions  of 
gravity  anomalies  or  of  deflections  of  the  vertical?  A reasonable 
answer  to  these  questions  requires  some  statistical  theory  of 
the  anomalous  gravitational  field.  But  is  this  field  really  a 
stochastic  phenomenon?  Such  questions  seem  to  motivate  research 
into  the  statistical  foundations  of  collocation. 

Least-squares  collocation  has  its  roots  in  many  fields: 

1.  Least-squares  estimation; 

2.  Prediction  theory  of  stochastic  processes; 

3.  Approximation  theory; 

4.  Functional  analysis,  especially  the  theory  of  Hilbert 
spaces  with  kernel  functions; 

5.  Potential  theory; 

6.  Inverse  and  improperly  posed  problems. 

All  of  these  "many  facets  of  collocation"  present  relevant 
aspects  which  must  be  taken  into  account  in  a complete  and 
balanced  treatment. 

The  relation  to  the  theory  of  inverse  problems  is  clear: 
our  data  are  functionally  related  to  the  gravitational  field;  to 
determine  this  field  from  the  data,  we  must  somehow  invert 
those  functional  relations.  Now  the  gravity  field  requires 
infinitely  many  parameters  for  its  full  determination;  the 
number  of  measurements,  however,  is  essentially  finite.  There- 
fore, we  have  an  improperly  posed  problem.  To  get  a unique 
solution,  we  must  impose  additional  conditions,  which  may  have 
the  form  of  a least-squares  principle  or  of  a norm  in  Hilbert 
space . 
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Historically,  collocation  has  developed  from  least- 
squares  prediction  of  gravity  anomalies,  which  is  an  application 
of  the  prediction  theory  of  stochastic  processes.  Hence, 
statistical  considerations  have  played  an  essential  role  in 
collocation  from  the  very  beginning. 

Also,  the  relation  to  classical  least-squares  adjust- 
ment has  soon  been  noted . In  fact,  collocation  models  bear  formal 
resemblance  to  conventional  adjustment  models.  The  characteristic 
difference,  however,  is  the  infinite  number  of  parameters 
necessary  to  fully  characterize  the  gravitational  field.  This 
fact  furnishes  an  essential  link  to  stochastic  processes  and  to 
infinite-dimensional  Hilbert  spaces. 

Least-squares  estimation  and  stochastic  processes  give 
a very  convenient  mathematical  formalism  and  terminology.  They 
also  provide  the  basis  for  a statistical  interpretation  of  the 
results,  essential  for  feasibility  studies. 

The  practical  success  of  the  statistical  treatment  of 
collocation  has  sometimes  overshadowed  its  equally  significant 
analytical  aspects,  especially  the  fact  that  there  is  a clean 
analytical  structure  underlying  it.  This  mathematical  structure 
is  based  on  the  harmonic  character  of  the  anomalous  gravitational 
field  and  on  the  fact  that  all  quantities  of  this  field  can  be 
expressed  as  linear  functionals  of  the  anomalous  potential.  The 
analytical  character  of  collocation  is  best  brought  out  by 
approaching  it  from  the  standpoint  of  approximation  theory, 
working  in  a Hilbert  space  with  a kernel  function. 

These  two  aspects,  the  statistical  and  the  analytical 
aspect,  are  both  indispensable  and  mutually  complement  each 
other.  This  fact,  evident  already  in  the  fundamental  paper 
(Krarup,  1969),  seems  to  be  generally  agreed  upon,  although 
there  is  some  controversy  on  details,  as  may  be  seen  from  the 
papers  collected  in  (Moritz  and  Siinkel , 1978) ; cf.  also  (Dermanis, 
1976) . 


3 


A literal  interpretation  of  the  anomalous  gravitational 
field  as  a stochastic  process  has  encountered  two  objections. 
First,  there  is  only  one  Earth;  a probability  space  of  many 
possible  earths  is  logically  unobjectabl e , but  appears  unnatural, 
since  all  realizations  except  one  (the  real  Earth)  are  un- 
observable. Secondly,  Lauritzen  (1973)  has  proved  that  there  is 
no  ergodic  Gaussian  process,  harmonic  outside  a sphere.  This  has 
sometimes  been  misinterpreted  as  a proof  that  no  ergodic  process 
modelling  the  anomalous  gravity  field  exists  at  all,  so  that 
the  covariance  function,  in  principle,  cannot  be  estimated  from 
the  data.  In  fact,  however,  the  Gaussian  structure  enters 
essentially  into  Lauritzen's  proof,  ar.d  there  do  exist  non- 
Gaussian  ergodic  processes  suitable  for  collocation. 

In  the  present  report  we  shall  attempt  an  elementary 
discussion  of  possible  stochastic  processes  on  the  sphere  which 
are  suited  as  statistical  models  for  the  earth's  gravitational 
field.  As  a preparation,  we  shall  first  consider  stochastic 
processes  on  the  circle,  which  are  simpler  and  already  show 
essential  theoretical  features. 

We  shall  present  two  different  ergodic  stochastic  process 
models.  One  is,  in  a way,  a non-Gaussian  ergodic  analogue  of 
Lauritzen's  model;  there  is  an  underlying  probability  space  of 
infinitely  many  different  "sample  earths".  For  the  second  model, 
the  probability  space  is  rotation  group  space;  all  realizations 
differ  only  by  a rotation,  so  that  there  is,  in  fact,  only  one 
Earth.  This  model  is  extremely  simple:  it  has  been  called 
"trivially  ergodic"  in  (Moritz,  1973,  p . 70 ) . At  the  same  time, 
it  expresses,  in  a natural  way,  the  homogeneity  and  isotropy  of 
the  anomalous  gravitation  field,  which  is  usually  presupposed 
in  collocation.  This  second  model  allows  a formal  statistical 
analysis  (covariances  and  distributions)  of  the  terrestrial 
gravitational  field  even  if  we  reject  the  interpretation  of  this 
field  as  a stochastic  phenomenon  in  a genuinely  physical  sense. 
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1.  Stochastic  Processes  on  the  Circle 

The  anomalous  gravitational  potential  of  the  earth  is 
a harmonic  function  outside  the  earth's  surface.  Outside  a 
certain  sphere,  such  a function  is  uniquely  determined  by  its 
values  on  the  spherical  surface:  from  these  values  it  is  ob- 
tained by  solving  an  exterior  Dirichlet  problem.  To  any  con- 
tinuous function  on  the  sphere,  a harmonic  function  in  outside 
space  can  be  made  to  correspond  in  this  way.  Instead  of  studying 
the  behavior  of  a spatial  harmonic  function,  we  may  thus  in- 
vestigate the  behavior  of  a (rather  arbitrary)  surface  function 
on  a sphere. 

It  is  in  this  sense  that  the  earth's  external  anomalous 
potential  has  frequently  been  mathematically  described  by  a 
stochastic  process  on  a sphere.  The  earth's  surface  is  very 
nearly  a sphere.  Therefore,  homogeneity  and  isotropy  on  the 
earth's  surface  may  approximately  be  formulated  in  terms  of  the 
rotation  group.  This  also  accounts  for  the  usefulness  of 
spherical  harmonics,  which  are,  in  a natural  way,  related  to 
the  rotation  group. 

For  the  present  purpose,  almost  all  essential  features 
are  preserved  if  we  consider  functions  in  the  plane  that  are 
harmonic  outside  a circle  instead  of  functions  in  space  that 
are  harmonic  outside  a sphere,  and  stochastic  processes  on  the 
circle  instead  of  stochastic  processes  on  the  sphere.  Further- 
more, this  reduction  of  dimensionality  essentially  simplifies 
the  problem  and  makes  its  mathematical  structure  easy  to  under- 
stand. Therefore,  we  shall  start  with  the  study  of  stochastic 
processes  on  the  circle. 

A continuous  and  continuously  differentiable  function 
f(t)  on  the  unit  circle  0 < t < 2*  can  be  expanded  into  a / 

uniformly  convergent  Fourier  series  (Smirnow,  v.II.  p . 4 1 7 ) : 


oo 


f(t)  s I (akcoskt  + bksinkt) 
k=o 
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where  ak  and  bK  are  coefficients;  since  sinkt  = 0 for 
k = 0 , bQ  is  arbitrary  and  will  be  put  equal  to  zero. 

This  representation  defines  f(t)  also  for  arbitrary 
real  t (-  » < t < »)  as  a periodic  function: 


f(t  + 2 k tt ) = f(t)  , k = 1,  2,  3, 


(1-2) 


; In  view  of  the  well-known  orthogonality  relations  of 

the  trigonometric  functions: 


/ coskt  costt  dt  = 0 if  k t i , 
o 

2 IT 

/ sinkt  sins-t  dt  = 0 if  k t i , 
o 

2 IT 

/ coskt  sinjtt  dt  = 0 always  , 
o 

2ir  „ 2 it  - 

/ cos  kt  dt  = /sin  kt  dt  = tt  if  k > 0 , 


(1-3) 


(1-4) 


the  coefficients  of  the  series  (1-1)  are  given  by 


i 2tf 

7-  Jf ( t)dt  , 


. ZTT 

= 7 Jf(t)coskt  dt  if  k > 0 , 

1 2ir 

= - /f(t)sinkt  dt  if  k > 0 . 


The  function,  defined  in  the  xy-plane  outside  and  on 
the  unit  circle. 


«■ 
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f(x,y)  = I r~k(a  coskt  + b sinkt)  (1-5) 

k=0  k k 

wi  th 

r = 1/x2  + y2  , t = arctan  ^ (1-6) 

being  polar  coordinates,  reduces  on  the  unit  circle  r = 1 to 
(1-1)  and  is  readily  seen  to  be  harmonic  for  r > 1 , satis- 
fying Laplace's  equation 

^ * *4  * 0 . (1-7) 

ax'*  ay* 

We  thus  have  a very  simple  one-to-one  relation  between  the 
function  (1-1)  defined  on  the  unit  circle  and  the  harmonic  func- 
tion (1-5)  defined  outside;  it  will  therefore  be  sufficient  in 
the  sequel  to  limit  our  study  to  (1-1). 

A stochastic  process,  or  random  function,  on  the  circle 
is  a function  f ( t ,u) ) which  depends,  in  addition  to  t , on 
a parameter  a>  which  represents  a "random  choice".  For  any 
fixed  value  u = u1  we  get  a function  f(t,u>1)  of  t only, 
which  under  the  above-mentioned  assumptions  has  the  form  (1-1); 
different  give  different  functions  of  t of  form  (1-1), 

which  are  considered  as  different  "realizations"  of  the  random 
process  f(t,u>)  . 

For  instance,  w may  denote  the  numbers  1,  2,  3,  4,  5, 

6,  so  that  f(t,oi)  denotes  6 functions  of  form  (1-1).  By 
throwing  a die  we  can  determine  w (e.g.,  u1  = 5)  and  the 
function  f ( t ,oj j ) associated  with  it;  this  will  explain  the 
term,  random  function. 

More  generally,  u>  is  a point  in  some  probability  space, 
or  sample  space,  n . In  this  space  we  define  a measure,  such 
that  measurable  subsets  of  n are  associated  with  events,  the 


X, 
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measure  of  a subset  denoting  the  probability  of  the  corresponding 
event.  The  measure  of  ft  itself  is  1. 

Let  us  illustrate  this  well-known  fact,  which  can  be 
found  in  any  textbook  on  probability  (my  favorite  is  (Feller, 
1957,  1966))  by  means  of  the  example  just  given,  the  throw  of 
a die.  Probability  space  ft  is  the  set  of  the  six  integers 
(1,  2,  3,  A. , 5,  6 } . Any  of  these  integers,  say  4,  forms  a 
subset  of  ft  , denoted  by  {4}  . This  subset  corresponds  to  the 
event  of  throwing  the  face  "4".  To  each  of  the  subsets  (1}, 

(2),  ...»  {6}  we  associate  the  same  measure  1/6.  The  event 
of  throwing  a "2"  or  a "4"  corresponds  to  the  sum  of  the  sets 
{2}  and  {4}  and  has  probability  equal  to  the  sum  of  the 
individual  probabilities: 


The  event  of  throwing  a "1"  or  a "2"  or  a "3"  or  a "4"  or  a "5" 
or  a "6"  has  the  probability 


1 

F 


. 1 . 1 + 
+ F + F + 


1 

F 
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that  is,  certainly,  as  it  must  be  from  an  intuitive  point  of 
view:  it  is  certain  that  one  of  the  faces  from  "1"  to  "6"  will 
show  up.  This  illustrates  the  intuitive  reason  for  demanding 
that  the  total  measure  of  ft  is  1. 

In  this  simple  example  we  have  6 possible  choices,  or 
"sample  points".  In  more  relevant  case  we  need  infinitely  many 
possible  choices,  corresponding  to  a more  general  probability 
space  ft  . 

Let  us  return  to  our  case  of  a random  function  on  the 

circle 


f = f ( t ,U)  ) , 0 < t < 2 TT 

w 6 ft  > 


d-8) 
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ft  denoting  a general  probability  space,  which  will  be  special- 
ized later  on.  To  get  these  simple  but  basic  concepts  firmly 
fixed  in  our  mind,  let  us  state  again  the  meaning  of  the  two 
arguments  t and  u , using  slightly  different  terms. 

The  variable  t is  the  “space  variable",  defining 
position  in  actual  physical  space.  This  becomes  immediately 
evident  on  taking  into  account  that  the  circle  is  a simplified 
analogue  to  the  terrestrial  sphere,  so  that  a point  on  the 
circle,  defined  by  t , corresponds  to  a point  on  the  earth's 
surface . 

On  the  other  hand,  u , so  to  speak,  describes  chance: 
it  defines  a random  choice.  In  statistical  mechanics,  the 
probability  space  ft  is  called  phase  space;  we  shall  sometimes 
find  this  terminology  convenient  and  call  u a phase  variable. 
Anyway,  w serves  a kind  of  "random  label"  to  distinguish  one 
realization  (or  sample  function)  f(t,a>1)  of  our  stochastic 
process  from  another  realization  f(t,oi2)  , both  sample  func- 
tions being  functions  of  t only,  since  w 1 or  «2  are  con- 
stants . 

Generally  speaking,  a quantity  depending  on  u>  is 
called  a random  variable.  This  explains  the  name,  random  func- 
tion, for  a function  f(t,u>)  of  t that  depends,  in  addition, 
on  "chance"  u . 

Let  us  expand  such  a random  function  on  the  circle  into 
a Fourier  series  (1-1)  with  respect  to  t . We  have 

oo 

- i r*t(»  )coskt  + b (w)sinktl  ; (1-9) 

k=0  K K 

clearly,  the  coefficients  ak  and  bR  will  now  be  random, 
variables  depending  on  u . By  (1-4)  they  are  given  by 


(1-1 


= 77  /^(t,uj)dt  , 

o 

i 2ir 

a.  (u>)  = T /f (t,u)coskt  dt  , k > 0 , 

K * o 

and  similarly  for  bk(u>)  • 

2 . The  Covariance  Function 

Consider  the  values  of  a random  function 
different  positions,  t and  t + s (Fig. 1 ) and 
product: 


f ( t ) f ( t+s ) ; 

t-0 


Figure  1.  Positions  on  the  circle 
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f at  two 
form  their 

(2-1 


l 


i 


the  dependence  on  u>  will  always  be  understood  even  if  not 
explicitly  written.  A suitably  defined  average  of  the  product 
(2-1)  is  nothing  else  than  the  covariance  function  corresponding 
to  the  random  function  f(t)  = f(t,u>)  ; it  depends  on  the  dis- 
tance s and,  possibly,  also  on  t and  <o  . 

For  random  functions,  the  natural  definition  of  the 
average  is  in  terms  of  the  statistical  expectation  E : 

C(s,t)  = E( f ( t) f ( t+s ) } 

= / f(t,u)f(t+s,u)dn  ; (2-2) 

ft 

E is  defined  as  an  integral  over  probability  space  ft  . This 
definition  presupposes  that  the  random  function  itself  has  zero 
expectation: 

E{ f } = /f (t,w)dft  = 0 . (2-3) 

ft 

By  (1-10)  this  implies 

E{ a } = 0 = E{ b } for  all  k . (2-4) 

JC  K 

We  substitute  the  Fourier  expansion  (1-9)  into  (2-2) 

and  get 

oo 

C(s,t)  = E{  l [akcoskt  + bksinkt]  • 
k=0 

OO 

• l [aAcost(t+s)  + b^ sint (t+s)]  } 

1 = 0 

* E{ n [ava» coskt  cost(t+s)  + 

kt 


+ bj^b^sinkt  siru(t+s)  + 
+ akb4coskt  sint(t+s)  + 
+ bj^a^sinkt  cos*(t+s)]} 


(2-5) 


The  formal  multiplication  of  the  two  Fourier  series  is  justi- 
fied since,  by  our  assumption,  they  are  uniformly  convergent. 
For  the  same  reason,  we  can  perform  the  integration  E term 
by  term. 

We  shall  now  make  the  fundamental  assumption  that  the 
Fourier  coefficients  are  all  statistically  uncorrel ated , that 
is,  that  all  covariances  between  different  coefficients  vanish; 


E{\V 

= 0 

if 

k i t , 

i E<W 

= 0 

if 

k i i , 

(2-6) 

II 

E(akbi} 

= 0 

always  . 

We  further  assume  that  the  variances  of  a,  and  b,  , for  each 

k k 

k , are  equal : 


E(a^)  - E<b*, 


Then  (2-5)  becomes 


(2-7) 


C(s,t)  = l [E{a^}coskt  cosk(t+s)  + 
k=0 

Efb^lsinkt  sink(t+s)]  . 


(2-8) 


•wiiw^T 
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In  view  of  (2-7)  and  of  the  identity 

coskt  cosk(t+s)  + sinkt  sink(t+s)  = cosks 
this  finally  reduces  to 


C(s)  = l ckcosks  , 

k=0 


(2-9) 


which  shows  that  the  covariance  function  then  depends  on  the 
distance  s only. 

The  Empirical  Covariance  Function.-  In  practice,  one 
frequently  has  only  one  realization  of  a stochastic  process 

f(t)  =•  f(t,w)  , u = const.  (2-10) 

The  question  is  whether  it  is  possible  to  estimate  the  co- 
variance  function  using  this  one  sample  function  only. 

In  this  case  we  cannot  form  the  statistical  expectation 
E , the  "phase  average"  (if  probability  space  si  is  denoted 
as  phase  space);  instead,  we  form  an  average  over  t , the 
"space  average"  M (for  stochastic  processes  on  the  real  line, 
- » < t < « , t may  be  interpreted  as  time,  so  that  M will 
be  a "time  average").  The  space  average  M of  (2-1)  is  defined 
as 


r (s)  = M{f (t)f (t+s)}  = ^ /f(t)f(t+s)dt  , (2-11) 

o 


f(t)  being,  as  always,  understood  as  a periodic  function  (1-2). 
The  function  r(t)  is  called  the  empirical  covariance  function. 


In  analogy  to  (2-3)  we  have  the  condition 


M{  f ( t ) } = / f ( t ) d t = 0 , 


which  means  by  (1-4)  that 


a = 0 . 
o 


(2-12) 


(2-13) 


This  is  not  an  essential  restriction  since  we  can  always  re- 
place f(t)  by  f(t)  - aQ  , for  which  the  zero-order  coefficient 
is,  in  fact,  zero.  Thus  we  may  assume  (2-13)  to  hold. 

Then  the  sum  in  (1-1)  begins  with  k * 1 , and  sub- 
stituting this  series  into  (2-11)  we  get 


2n  “ 


1 -c.  II  — 

r(s ) s Yn  f [akcoskt  + bksinkt] 


o k=_ 


• l fa  cos£(t+s)  + b sin«(t+s)l  dt 
£-1“  * 1 


oo  00  2lt 


1 ~ 4.  M 

77  l X /[a^coskt  cost(t+s)  + 


k= 1 £=1  O 


+ bkbAsinkt  sinfc(t+s)  + 

+ akb£coskt  sin£(t+s)  + 

+ bRaJlsinkt  cost(t+s)]dt  , 

the  formal  operations  (series  multiplication  and  termwise  inte- 
gration) are  again  justified  by  uniform  convergence. 

The  orthogonality  relations  (1-3)  give  at  once: 


1 00  o2lr 

r(s)  = 77  I [ak  /coskt  cosk(t+s)dt  + 


+ /sinkt  sink(t+s)dt  + 


+ a.b  /coskt  sink(t+s)dt  + 


(2-14) 


+ akbk  /sinkt  cosk(t+s)dt~j  , 


since  all  products  of  trigonometric  functions  for  k / i 
vanish  after  integration. 

We  further  have 


/ coskt  cosk(t+s)dt 
o 

2ir  2 

= /(cos  kt  cosks 
o 

* it  cosks  , 


- coskt  sinkt  sinks)dt 


again  by  (1-3),  and  similarly 


/sinkt  sink(t+s)dt  = w cosks 


I 


Final ly. 


!lt  _ 

/[coskt  sink(t+s)  + sinkt  cosk(t+$)Jdt 


■ /sink(2t+s)dt  * 0 . 
o 


: 
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Hence  (2-14)  reduces  to 


oo 

r ( s ) * y l (ak  + bk)cosks  . 


(2-15) 


This  is  the  Fourier  expansion  of  the  empirical  co- 
variance  function.  In  view  of  (2-10),  the  function  f ( t ) and 
its  Fourier  coefficients  ak  and  bk  depend  on  o>  : 


ak.  = * bk  = bk (<*> ) . 


(2-16) 


Hence,  also  r(s)  depends  on  u , so  that  (2-15)  can  be  written 
more  explicitly: 


r(s,o>)  = l y (w)cosks  , 


where 


(2-17) 


Yk(“)  = C®k  ^ + (“>)[]  * 


(2-18) 


Let  us  now  compare  the  empirical  covariance  function 
(2-15)  or  (2-17)  with  the  true  covariance  function  (2-9). 
Forming  the  expectation  E of  (2-18)  we  have 


E{yk>  = E{ Yk (w ) } = y E(ak>  + ^ E{bk>  » 


III 


so  that  by  (2-7) 


E{V  " ck  * 


(2-19) 


tMiiNS 


The  expectation  of  (2-17)  is 


so  that 


(2-20) 


the  expectation  of  the  empirical  covariance  function  is  the  true 


covariance  function.  In  statistical  terms,  the  empirical  co- 
variance  function  is  an  unbiased  estimate  of  the  true  covariance 
function . 

It  would  be  particularly  desirable  if  the  empirical  co- 
variance  function  is  identical  to  the  true  covariance  function, 
or  if  the  two  function  are  equal  at  least  for  almost  all  u> 

(that  is,  for  all  u>  with  the  possible  exception  of  a set  of 
measure  zero).  In  this  case,  the  covariance  function  can  be 
exactly  estimated  from  one  realization  of  the  stochastic 
process  f(t,w)  , that  is,  from  one  sample  function  u = const. 
This  is  the  case  of  ergodicity. 

This  name  has  been  taken  from  statistical  mechanics, 
where  it  means  that  a time  average  is  the  same  as  the  corre- 
sponding phase  average.  In  our  case,  the  space  average  M of 
f(t)f(t+s)  should  be  equal  to  the  phase  average  E of  this 
product. 

Obviously,  ergodicity  is  a very  special  case,  and  the 
question  arises  whether  it  is  possible  at  all.  This  question 
will  be  answered  positively  in  the  next  section. 


3.  Ergodic  Processes  on  the  Circle 


The  case  in  which  the  empi rical  covariance  function 
r(s,«)  coincides,  for  almost  all  u , with  the  true  one, 

C(s)  , has  been  called  ergodicity  in  the  preceding  section.  By 
comparing  the  coefficients  of  the  respective  Fourier  expansions 
(2-9)  and  (2-15)  we  get  the  necessary  and  sufficient  condition 
for  ergodicity: 

ak(“)  + = 2ck  (3-1) 

for  almost  all  w . 

The  meaning  of  this  condition  should  be  carefully  kept 
in  mind.  The  coefficients  ck  , defined  by  (2-6),  are  given 
nonrandom  constants.  The  coefficients  a.  and  b.  on  the  left- 

k k 

hand  side  are,  however,  functions  of  w and  hence  random 

variables.  Thus  the  condition  (3-1)  is  certainly  very  restrictive. 

It  should  be  recalled  that  we  have  derived  (3-1)  under 
the  assumption  of  uniform  convergence  of  the  Fourier  series 
for  f(t,w)  . This  assumption  is  not  essential;  for  a proof 
under  more  general  conditions  ( integrabi 1 i ty)  see  (Zygmund, 

1968,  pp. 36-37). 

Lauritzen's  Theorem. - In  particular,  it  is  impossible 
to  satisfy  the  ergodicity  condition  by  a stochastic  process 
defined  by  (1-9)  with  ak(u>)  and  bk(w)  being  uncorrelated 
and  normally  distributed  (Gaussian)  stochastic  variables  of  zero 
expectation.  This  has  been  proved  by  Lauritzen  (1973,  p.65) 
by  explicitly  calculating  the  variance  of  the  empirical  co- 
variance  function  r(t,u)  and  showing  that  it  is  non-zero. 

(For  ergodic  processes  this  variance  is  evidently  zero.) 

For  us,  Lauritzen's  theorem  is  an  obvious,  almost 
elementary  consequence  of  (3-1).  For  Gaussian  random  variables, 
uncorrel atedness  is  equivalent  to  statistical  independence. 


18 


Hence  (2-6)  implies  that  all 


are  statistically 


independent  random  variables.  If  the  functions  a (<*> ) and 
bk(uj)  can  vary  independently  of  each  other,  then  (3-1)  can 
be  violated  at  will.  Eq.  (3-1)  would  only  be  satisfied  if 


ak(w)  = const.,  bk(w)  = const. 


(3-2) 


m 


for  almost  all  w , which  is  incompatible  with  zero  expectation 
(2-4).  These  contradictions  prove  the  theorem. 

Lauritzen's  theorem  may  be  concisely,  though  somewhat 
loosely,  formulated  thus:  a Gaussian  random  process  on  the 
circle  cannot  be  ergodic.  Looking  for  an  ergodic  process,  we 
must,  therefore,  consider  non-Gaussian  processes.  The  ak  and 
bk  will  be  uncorrelated,  but  not  necessarily  statistically 
independent.  It  is  known  that  statistical  independence  implies 
uncorrel atedness ; the  converse  is  true  only  for  normal  processes. 

Ergodic  Processes:  First  Example.-  Let  the  coefficients 


and  bk  , for  different  k , be  statistically  independent; 
for  the  same  k , ak  and  bk  will  only  be  uncorrelated,  in 


agreement  with  the  third  equation  of  (2-6).  To  satisfy  the 
ergodicity  condition  (3-1  * take 


a,  = /2c,  co 


b,  = /2c.  sino. 
k k k 


(3-3) 


where  u is  a random  variable  uniformly  distributed  in  the 
k 

interval  0 < u>k  < 2*  . This  means  that  the  probability  density 
$(a>  ) of  u has  the  form  of  Fig. 2.  Geometrically,  a and 

K JC  Jv 

bk  are  represented  in  Fig. 3,  We  have  a random  vector  with 
components  (a  ,b  ) of  fixed  length  /2c  but  with  randomly 
variable  azimuth  <o  . The  end  point  of  this  vector  thus  des- 

k 

cribes  a circle.  Any  point  of  the  circle  corresponds  to  a random 
choice  of  wk  . The  probability  that  the  end  point  of  the  vector 
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\Y 


falls  onto  the  arc  AB  is  proportional  to  the  arc  length, 
namely  (b-ci)/2it  , corresponding  to  the  shaded  area  in  Fig. 2; 
this  is  the  meaning  of  uniform  distribution.  Obviously,  the 
probability  that  the  end  point  lies  somewhere  on  the  circle  is 
2tt/ 2tt  = 1 , namely  certainty,  as  it  should  be. 

Thus  the  coefficients  aR  and  bk  are  clearly  statisti 
cally  dependent,  but  are  they  still  uncorrelated?  We  have 


E{akbk}  = K K>bkK^kK>d“k 

o 


(3-4) 


where 


<f>k(“>)  = 


77  * 0 < o)k  < 2*  , 


(3-5) 


so  that  by  (3-3) 


E(akV  = /2ckcos»ks1ri„lt  ^ d»k 


C,  2 TT 
k t 


7Z  ^Sin2“kda)k  = 0 ; 


(3-6) 


hence  a and  b are,  in  fact,  uncorrelated.  This  provides  a 

Jv  K 

simple  geometrical  illustration  of  how  two  random  variables  can 
be  uncorrelated  without  being  statistically  independent. 

For  different  k , the  vectors  (ak»bk)  have  been 
supposed  to  be  independent.  This  means  that  the  a>k  are  uni- 
formly distributed,  independent  random  variables.  Consider  two 
co k , say,  u>2  and  a>3  . Each  u>k  varies  from  0 to  2 it  , 
that  is,  over  the  (unit)  circle.  Since  the  joint  probability 
space  of  two  independent  random  variables  is  the  Cartesian 
product  of  the  two  individual  probability  spaces,  the  joint 
probability  space  of  u>2  and  a>3  is  the  Cartesian  product  of 

two  circles.  The  joint  probability  space  of  (a^,  u>2 wn) 

is  the  product  of  n circles. 
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The  Fourier  series  of  the  random  function  f(t,u>)  in- 
volves infinitely  many  coefficients  ak  and  t>k  . The  probabi- 
listic event  of  sorting  out  one  sample  function  thus  requires 
infinitely  many  independent  choices  of  u^,  a>2,  u,3,...  The 
probability  space  for  f(t,oj)  is,  therefore,  the  Cartesian 
product  of  infinitely  many  circles,  or  u>  represents  the  in- 
finite vector 


= (<*>!  . 


<1>2*  “3  * • • • ) » 


(3-7) 


each 


wk  being  independently  uniformly  distributed. 


Finally  we  prove  that  if  one  sample  function  of  our 
present  ergodic  process  has  a uniformly  convergent  Fourier 
series,  then  this  will  hold  for  ail  sample  functions  of  this 
process.  Since  the  absolute  values  of  sine  and  cosine  are,  at 
most,  equal  to  1,  the  Fourier  series  (1-1),  with  aQ  = 0 , has 
the  majorant 


l (laJ  + IbJ) 


(3-8) 


Convergence  of  this  majorant  series  is  clearly  sufficient  for 
uniform  convergence  of  our  Fourier  series;  that  it  is  also 
necessary  is  a consequence  of  the  Theorem  of  Denjoy-Lusin 
(Zygmund,  1968,  p.232). 

Since 


/a2  + bS  < | a | + |b|  < 2/?  TV  , 


(3-9) 


convergence  of  (3-8)  is  logically  equivalent  to  the  convergence 
of 


£yak + bk  “,s/frk  • 


(3-10) 
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which  is,  therefore,  also  a necessary  and  sufficient  condition 
for  the  uniform  convergence  of  our  Fourier  series.  Therefore, 
uniform  convergence  of  the  Fourier  series  of  one  sample  function 
implies  convergence  of  the  right-hand  side  of  (3-10).  Since  this 
right-hand  side  does  not  depend  on  w , the  left-hand  side  of 
this  equation  must  converge  for  all  a>  , which  implies  uniform 
convergence  of  the  Fourier  series  of  the  sample  functions  for 
any  oi  , which  was  to  be  shown. 

The  uniform  distribution  on  a circle  is  even  simpler 
than  a normal  distribution.  Furthermore,  the  "probability  circle" 
0 < a)  < 2tt  seems,  somehow,  to  be  a natural  counterpart  of  the 
"space  circle"  0 < t < 2ir  . Thus  the  present  simple  example 
seems  to  be  a quite  natural  model  for  a stochastic  process  on 
the  circle,  more  natural  than  any  Gaussian  model;  furthermore 
it  is  ergodic.  The  next  example  is  still  simpler. 

Ergodic  Processes;  Second  Example.-  We  now  take  u it- 
self as  a random  variable  uniformly  distributed  in  the  interval 


0 < 


2 TT 


(3-11) 


or,  what  is  the  same,  on  the  unit  circle.  Thus,  in  the  random 
function  f(t,u>)  , both  variables  t and  u>  now  range  over  a 
unit  circle,  the  circle  for  t representing  "ordinary  space" 
and  the  circle  for  u representing  "probability  space". 

We  now  take 


f(t,u>)  = f(t+w)  . 


(3-12) 


Let 


f(t,a>)  = f ( t+u> ) represents  simply  a rotation  of  the  circle,  or 
of  the  function  f(t)  , by  the  angle  u>  (Fig. 4). 


\f(t*w) 


Figure  4.  The  rotation  f(t)  *=>  f(t+u>) 


We  may  also  write 


f ( t+oi ) * Rwf(t)  , 


(3-14) 


where  the  operator  means  rotation  by  the  angle  u>  . In 

other  terms,  we  may  identify  our  probability  space  (3-12)  with 
rotation  group  space.  In  fact,  in  the  plane,  the  rotation  group 
is  one-dimensional,  being  characterized  by  one  angle  u>  . 

The  functions  f(t,a))  differ  from  each  other  only  by 
a rotation;  they  are  not  essentially  different  (Fig. 4).  This 


model,  therefore,  is  suited  to  represent  the  case  in  which  there 
is  only  one  realization  f(t)  and  we  wish  to  use  the  mathe- 
matical techniques  of  stochastic  processes;  this  is  justified 
in  the  case  of  homogeneity,  if  a rotation  by  u>  gives  a 
physically  equally  meaningful  situation,  so  that  f(t+w)  in- 
stead of  f(t)  would  be  physically  equally  possible. 

We  assume  (2-12)  to  hold  for  our  initial  f(t)  . Then 
(2-3)  becomes 


, 2 IT 

E{ f ( t ,<u ) } = £7  Jf(t,«)d« 


(3-15) 


/ f ( t+u) ) dw 


We  change  the  integration  variable  by  putting  ( t is  constant 
with  respect  to  integration1.) 


t + u>  = u , du  = du  , 


(3-16) 


obtaining 


t + 2 IT 


E{f(t,u>)}  = / f (u)du 

= 77  /f(u)du  = 0 


(3-17) 


by  (2-12),  so  that  (2-3)  is  satisfied. 

Now  the  covariance  function  (2-2)  becomes 


i 2* 

C(s,t)  = -k — J f ( t+co ) f ( t+ S+ui ) du>  . 


(3-18) 


The  substitution  (3-16)  transforms  this  integral  into 
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Ar  /f(u)f (u+s)du 

o 


- 77  /f ( t ) f ( t+s ) dt 

= r(s) 


(3-19) 


by  (2-11).  Therefore,  in  this  model,  the  empirical  covariance 
function  coincides  with  the  true  covariance  function;  the 
process  is  ergodic.  In  fact,  the  "phase  average"  E is  seen 
to  coincide  with  the  "space  average"  M . Since  E can  be 
transformed  into  M by  a simple  change  of  variables,  the 
process  under  consideration  is  trivial ly  ergodic. 

This  is  obviously  a very  simple  situation,  but  an 
important  one.  We  shall,  therefore,  try  to  understand  it  better 
by  studying  the  spectral  representation,  that  is,  the  Fourier 
series  expansion. 

We  shall  denote  the  Fourier  coefficients  of  the  initial 


representation  f(t)  = f(t,0)  by  ak  and  bk  . The  coeffi- 
cients ak(io)  and  bk(o>)  of  f(t,a>)  are  then  given  by  (1-10) 


We  have 


, 2ir  1 2tt 

(oj)  = -x—  /f(t+u)dt  = ■*—  /f(t)dt  = 0 


(3-20) 


by  (2-12).  For  k > 0 we  get 


, 2n 

av(“)  = ~ /f(t+w)coskt  dt 


(3-21) 


and  substituting  ( u>  is  constant  with  respect  to  integration!) 


t + u = v , dt  = dv  , 


(3-22) 
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we  have 


a.  (u>)  = T Jf (v)cosk(v-u)dv 
K o 

, 2n 

= 1 I f(v)(coskv  coskw  + sinkv  sinkw)dv 
* o 

= coska)  • - /f(v)coskv  dv  + 


(3-23) 


- ZTl 

+ sinka)  • - /f (v) sinkv  dv 


or,  by  (1-4), 


ak(w)  = akcoskt*>  + bksinkw  . 

In  exactly  the  same  way,  replacing  coskt  by  sinkt 
we  get 

b.  (<*) ) = - a.sinko)  + bkcoska>  . 


(3-24) 


sink(v-w) 


(3-25) 


Let  us  now  evaluate  (2-6),  using  (3-24)  and  (3-25).  We 


1 2ir 

E{ak(w)aJl(a))}  = -2J  Iak  (“ ) ai  (“ ) 


1 Z IT 

= J(akcoska>  + bksinkw)  * 


(a^cosno)  + b^sintuOduj  . 


On  termwise  multiplication  and  integration,  using  the  ortho- 
gonality relations  (1-3),  we  readily  obtain  the  value  zero  if 
k i i . Proceeding  similarly,  we  see  that  all  orthogonality 


! 


F 


relations  (2-6)  are  satisfied. 
We  further  obtain 


EUk(u>n 


2 IT 


1 2 
= 7-  /(acosku  + b sinku)  du 

o K 

* - bk>  • 


(3-26) 


In  the  same  way. 


E(b^(„))  = ♦ b*>  , 


(3-27) 


which  is  independent  of  w , so  that  (2-7)  is  satisfied  with 


ck  * 7<ak  + bk>  • 


(3-28) 


We  finally  compute,  using  (3-24)  and  (3-25), 


a (u)2  + b (u)2  = (a  cosku  + b sinku)2  + 


+ (-  a sinku  + b cosku) 

Jv  JC 


a2  a.  h2 
9k  + bk  • 


Thus 


ak(u)2  + bk(u)2  = 2ck  , 

which  shows  that  the  ergodicity  condition  (3-1)  is,  in  fact, 
satisfied . 


Let  us  finally  compare  this  model  with  our  first  ergodic 
model.  In  the  first  model,  probability  space  is  the  Cartesian 

product  of  infinitely  many  circles  0 < a>k  < 2ir  (k  = 1,  2,  3,...) 

<d  being  the  infinite  vector  (3-7),  consisting  of  independent, 
uniformly  distributed  random  variables.  In  the  present  model, 
probability  space  is  simply  one  circle  0 < u>  < 2n  , u>  being 
a uniformly  distributed  one-dimensional  random  variable.  There- 
fore, in  the  first  model,  ak  and  a^  for  k + i , depending 
on  different  independent  random  variables  a>k  and  , are 
statistically  independent.  The  same  holds  for  ak  and  b^  , 

and  for  bk  and  b^  . For  the  same  k , aR  and  bk  are 

dependent  though  uncorrelated.  On  the  other  hand,  in  the  present 
model,  all  Fourier  coefficients  depend  on  the  same  variable  w ; 
therefore,  all  are  statistical  dependent,  but,  as  we  have  seen, 
any  two  different  coefficients  are  uncorrelated,  as  a consequence 
of  the  orthogonality  relations  (1-3)  for  trigonometric  functions. 


4.  Stochastic  Processes  on  the  Sphere 

Notations . - Our  preceding  considerations  about  stochastic 
processes  on  the  circle  can  be  translated  almost  literally  to 
the  sphere.  Instead  of  the  Fourier  series  (1-1)  we  have  the 
spheri cal -harmonic  series 


°°  n 

f (e  ,x  ) = l l a R (e  ,x  ) + b S (e  , a ) 

' * ' L L nm  nm ' ’ ' nm  nm ' ’ ' 

n=0  m=0  >- 


(4-1) 


where 


R (e,x)  = p (cose)cosmx  , 
nm'  ' nm ' ' 


s„m(0»x)  s Pnm(cose)sinrnx  , 


(4-2) 
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P (cose)  being  the  (conventional)  Legendre  functions  (cf. 
Heiskanen  and  Moritz,  1967,  p.29);  n and  m are  called  degree 
and  order,  respectively. 

To  simplify  the  notation,  let  us  put 


Snnj(e,A)  = » m=  1 , 2 , . . . , n , 


(4-3) 


!: 


so  that  any  R with  negative  second  subscript  denotes  the 
nm 

corresponding  S , for  instance,  R_  = S_.  . Then  (4-1) 

nm  5 , -3  5 j 

may  be  simply  written  as 


■ l I - 

n=0  m=-n 


(4-4) 


if  for  the  coefficients  we  use  an  analogous  notational  con- 
vention: 


a — b , m 1,2,. . . , n 
n , -m  nm 


(4-5) 


Sometimes  it  will  be  useful  to  use  fully  normalized 
harmonics,  denoted  by  1?  and  S , or  by  IT  with 

nm  nm  nm 

- n < m < n , which  differ  from  the  conventional  harmonics  by 
a factor  and  are  normalized  by 


i-  //IT2  do  = 1 , 

4w  1 * nm 


(4-6) 


o denoting  the  unit  sphere  (ibid. , p . 3 1 ) . Spherical  harmonics 

are  orthogonal  functions:  if  we  integrate  the  product  of  any 

two  different  functions  R (or,  of  course,  IT  ) over  the 

nm  nm 

sphere,  we  get  zero.  The  fully  normalized  spherical  harmonics 
form  a system  of  orthonormal  functions: 
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W{  R R : 
nm  qp 


<5  6 

nq  mp 


where 

W * > = 37  //( * )d o 
a 


(4-7) 


(4-8) 


denotes  now  the  average  over  the  unit  sphere  and  is  the 

Kronecker  delta,  1 if  k = i and  0 otherwise. 

If  we  write  (4-4)  in  fully  normalized  harmonics. 


f(e.»>  - l l >„»,(«.>]  . 


n»0  m= -n 


then  the  coefficients  are  simply  given  by 


a-  = M{  fTT  } , 

nm  nm 


(4-9) 


(4-10) 


in  view  of  the  orthonormality;  these  equations  are  a shorthand 
notation  of  eqs.  (1-76),  ibid.  , p . 3 1 . 

As  a final  notational  convention  regarding  spherical 
harmonic  expansions,  we  introduce  the  two-dimensional  parameter 

t = (e , a ) (4-11) 

and  write  (4-9)  as 

f(t)  • ? l a T?  (t)  ; (4-12) 

' ' L nm  nm'  ' ' ' 

n=0  m=-n 

this  stresses  the  analogy  with  the  case  of  the  circle. 


The  stochastic  parameter  will  again  be  denoted  by 
os  € Q , a being  probability  space  with  total  measure  1.  Then 

f ( t ,u> ) 

will  denote  a stochastic  process  on  the  sphere.  The  expectation 
E is  agair  defined  by 

E{ • } = / ( - )dn  (4-13) 

fl 

as  an  average  over  probability  space,  or  phase  average,  as  opposed 
to  the  space  average  TT  defined  by  (4-8). 

In  analogy  to  (1-5),  there  is  a one-to-one  correspon- 
dence between  continuous  functions  on  the  sphere  and  harmonic 
functions  in  space:  the  spatial  function 


f (x.y.z) 


00 


= l 

n = 0 


nm 

y.n+  1 


R 


nm 


(e,x) 


(4-14) 


satisfies  Laplace's  equation  outside  a . Therefore,  there  is  a 
one-to-one  correspondence  between  harmonic  stochastic  processes 
in  space  and  stochastic  processes  on  the  sphere,  and  we  can 
limit  our  considerations  to  the  latter. 

Covariances ■ - We  again  assume  that  our  stochastic 
process  is  centered: 


E{ f ( t ,oi) } = 0 . (4-15) 

Then  the  covariance  C(t,u)  between  f(t,u))  and  f(u,aj)  at 
two  different  points  t and  u on  the  unit  sphere  a is,  as 
usual,  defined  by 


C(t,u)  = E{f (t)f(u)} 


t 


(4-16) 


the  dependence  on  u being  understood. 

As  in  the  circular  case,  we  shall  limit  ourselves  to 
continuously  differentiable  functions.  Then  the  spherical- 
harmonic  expansion  will  be  a uniformly  convergent  series 
(Kellogg,  1929,  p.259),  which  can  be  multiplied  and  termwise 
integrated . 

We,  therefore,  substitute  (4-12)  into  (4-16): 


C(t,u) 


E(  l l 

n=0  m=-n 


a 1?  (t) 

nm  nm ' ' 


multiply  and  integrate  termwise  with  respect  to  u>  (that  is, 
interchange  the  order  of  summation  and  integration),  obtaining 


C(t.u)  = l I l l E{a  a }K  ( t)TT  (u)  . (4-17) 

' ' l l l l nm  qp  nm ' ' qp ' ' ' ' 

n m q p ^ 

Let  us  now  assume  that  the  coefficients  a"  - J (w) 

nm  mu'  ' 

are  mutually  uncorrelated  random  variables: 


E{a  a } 
nm  qp 


(4-18) 


if  q / n or  p t m or  both,  and  that  E{a2  } is  the  same 

nm 

for  all  coefficients  of  degree  n , that  is,  for  all  m ; we 
put 


E{Inm>  = TnTT 


(4-19) 


Then  (4-17)  becomes 
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c<‘-u>  ■ l ZSTT  • (4-z(>) 

n=0  m=-n 

Now  we  make  use  of  the  well 'known  decomposition  formula  of 
spherical  harmonics  (cf.  Heiskanen  and  Moritz,  1967,  p.33,eq. 

(1-82')),  which  in  our  present  notation  takes  the  form 

Pn(cos<|>)  = 7~T  l Knm(t)tfnm(u)  , (4-21) 

m=-n 

with 

t = (0  ,X)  and  u = (6 ' ,X ' ) , (4-22) 

</>  being  the  spherical  distance  between  the  points  t and  u : 
cosiji  = cosecose'  + si  nesi  ne  ‘ cos  (x  ' -x ) , (4-23) 

i 

and  Pn(cos^)  denoting  the  (conventional)  Legendre  polynomial 
of  degree  n . Thus  (4-20)  reduces  to 

C(*)  ■ l CnPn(cos^,)  . (4-24) 

n=0 

Thus,  the  covariance  function  depends  only  on  the 
spherical  distance  ^ . This  is  the  important  case  of  homo- 
geneity and  isotropy;  it  is  seen  to  result  from  the  postulate 
that  the  variances  (4-19)  of  all  coefficients  a"  of  the 

nm 

same  degree  n are  equal . 

The  Empirical  Covariance  Function.-  If  there  is  only 
one  realization  of  the  stochastic  process,  we  cannot  directly 
compute  the  true  covariance  function  C defined  by  (4-16) 
and  expressed  by  (4-24).  We  may  again  try  to  compute  an  empirical 
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covariance  function  r by  replacing  the  phase  average  E by  a 
suitable  space  average  and  hope  that  r will  be  a good  estimate 
of  C ; if  possible,  r should  even  be  equal  to  C . 

In  view  of  the  homogeneity  and  isotropy,  we  must  inte- 
grate not  only  over  the  sphere  (homogeneity),  but  in  addition 
over  the  azimuth  (isotropy).  Therefore,  we  must  supplement  the 
average  M , defined  by  (4-8),  by  additionally  averaging  over 
the  azimuth  a . The  resulting  average  M may  be  defined  by 

. 2tt  tt  2tt 

M{  • > = Iff  ( • ) si nededxda  . (4-25) 

8tt  x = 0 6=0  ct=0 

The  geometric  situation  is  shown  by  Fig. 5. 


Figure  5 


Integration  over  rotation  group  space 


The  averaging  is  first  performed  over  the  circle  ^ = const.  , 
whose  center  t = (6, A)  is  then  made  to  vary  over  the  whole 
sphere  a . 

The  angles  A,  e,  a can  be  regarded  as  the  three 
Eulerian  angles  defining  a rotation  in  three-dimensional  space, 
that  is.  A,  e,  a are  the  coordinates  of  an  element  of  the 
rotation  group  or,  of  a "point"  in  "rotation  group  space". 
Therefore,  M will  be  called  a rotation  group  average. 

Hence  the  empirical  covariance  function  is  given  by 


r(*)  = M{f(t)f(u))  , 


(4-26) 


where  M is  defined  by  (4-25)  and  the  points  t and  u have 
the  spherical  distance  ip  , which  is  constant  with  respect  to 
the  integration. 

Because  of  the  way  in  which  the  average  M is  computed, 
the  empirical  covariance  function  will  depend  only  on  the 
distance  . It  can,  therefore,  be  expanded  into  a series  of 
Legendre  polynomials  of  <|<  : 


r(<J»)  = l y P (cosi|/) 


(4-27) 


Yn  can  be  expressed  in  terms  of  the  spheri cal -harmoni c 


coefficients  a of  the  same  n , in  full  analogy  to  (2-18). 
nm 

The  derivation  is  given  in  (Heiskanen  and  Moritz,  1967,  pp.257 
259);  the  result  is  eq.  (7-28),  1 oc . cit. , which  in  the  present 
notation  reads 


* i »; 


(4-28) 


Note  that  this  very  simple  expression  is  obtained  by  using  con 
ventional  harmonics  on  the  left-hand  side  and  fully  normalized 
harmonics  on  the  right-hand  side. 

Clearly,  Yn  » as  well  as  anm  , are  random  variables, 
that  is,  functions  of  to  , Their  expectation  is  given  by 

E{yn>  = E{Yn(u)}  = I E{a^m(w)}  . 


In  view  of  (4-19)  this  becomes 
E<*n>  " Cn  • 


(4-29) 


so  that 


E{r(*)}  = C(*)  , (4-30) 

exactly  as  in  the  circular  case  (2-20):  r(<>)  is  an  unbiased 
estimate  of  C(<|>)  . 


5 . Ergodic  Processes  on  the  Sphere 

For  an  ergodic  process,  coincides  with  C(^)  . 

The  ergodicity  condition,  correspond! ng  to  (3-1),  is 


l 


m* 


(to) 


= C 

n 


» 


(5-1) 


for  almost  all  <o  ; c is  Independent  of  to  . This  condition 
is  equivalent  to  Yn  = cn  . 


Lauritzen's  Theorem.-  Assume  that  ^nm(«)  are  normally 
distributed  (Gaussian)  random  variables.  For  Gaussian  variables, 
uncorrel atedness  is  equivalent  to  statistical  independence.  From 
our  basic  presupposition  (4-18)  it  thus  follows  that  the  anrn(u>) 
must  be  statistically  independent  of  each  other.  Then  the 
summands  on  the  right-hand  side  are  independent  functions  of  u>  , 
so  that  (5-1)  will  be  violated  for  almost  all  w . Loosely 
formulated:  a Gaussian  random  process  on  the  sphere  cannot  be 
ergodi c . 

The  present  simple  proof  of  Lauritzen's  theorem  suffers 
from  the  slight  logical  defect  that  (5-1)  has  been  derived  on 
the  assumption  that  our  stochastic  process  is  sufficiently 
smooth  (differentiable).  Since  Gaussian  random  variables  may 
take  arbitrarily  large  values,  the  convergence  of  the  correspon- 
ding spher i ca 1 -harmon i c series  cannot  be  guaranteed;  still  less 
are  we  sure  that  the  corresponding  realizations  will  all  be 
differentiable  (this  may  even  be  a practical  argument  against 
admitting  a Gaussian  process  as  mathematical  model  for  the 
terrestrial  gravity  field).  In  fact,  (5-1),  just  as  (3-1),  holds 
for  more  generally  assumptions,  but  we  have  not  proved  this 
because,  for  the  present  ergodic  models,  differentiability  can 
be  presupposed. 

Thus,  our  deduction  of  Lauritzen's  theorem  has  the 
character  of  a plausibility  argument  rather  than  of  a fully 
rigorous  mathematical  proof.  It  has,  however,  the  decisive 
advantage  of  showing  the  essential  statistical  situation  under- 
lying it,  and  the  fact  that  the  Gaussian  character  of  the  process 
is  essential  to  the  theorem:  only  for  Gaussian  distributions, 
uncorrel atedness  implies  statistical  independence.  This  simple 
fact  is  hidden  below  the  mathematical  intricacies  of  Lauritzen's 
(1973, p. 65)  proof  and  has  not  always  been  clearly  understood. 


We  shall  now  consider  two  examples  of  (non-Gaussian) 
stochastical  processes  on  the  sphere,  corresponding  to  the  two 
examples  for  the  circle  given  in  sec.  3. 

First  Example;  Uniformly  Distributed  Coefficients. - 
The  two  Fourier  coefficients  ak  and  bk  define  a two-dimen- 
sional vector  whose  end  point  lies  on  a circle  of  radius  /2cR 
(Fig. 3).  Similarly,  the  2n+l  coefficients  anm  ( n fixed, 

- n < m < n ) form  a (2n+l)-dimensional  vector 

a = fa  , a , ...  , a . , a„  1 (5-2) 

— L n,-n  * n,-n+l  n,n-l  n,n  J 

whose  end  point  lies  on  a sphere  of  radius  /c~  in  R2n+1 
(Euclidean  space  of  dimension  2n+l);  in  fact,  (5-1)  may  be 
written 

1 aj  2 ■ c„  . (5-3) 

Assume  now  that  different  realizations  of  the  stochastic 
process  correspond  to  different  positions  of  the  endpoint  of  a 
on  this  sphere.  In  other  terms,  if  e is  the  unit  vector 
corresponding  to  a , then 

a(w)  = /c  e(u>)  , ( 5-4 ) 

— ti- 

the unit  vector  being  a function  of  u>  : the  random  vector  a 
has  a random  direction  but  a constant  length,  in  complete 
analogy  to  (3-3).  The  random  directions  e(w)  are  uni  form! y 
distributed  in  our  (2n+l)-dimensional  space:  probability  is 
given  by  an  area  on  the  unit  sphere  in  this  space;  cf. (Feller, 
1967,  p.  68)  for  R3  . 

In  view  of  (5-3),  the  2n+l  coefficients  a of  the 
same  degree  n are  statistically  dependent,  but  they  are 
uncorrelated:  it  is  easy  to  see  that  (4-18)  holds  for  them.  In 
fact,  this  means  that,  for  two  different  components  of  the 
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vector  a or  of  the  vector  e , say  ei  and  e^  , the  integral, 
over  the  unit  sphere  a2n  in  R2n+i  * ^s  Pr°duct  eiej  is 
zero: 

/ eiejda2n  = 0 * ^5_5) 

°2n 

Denote  the  left-hand  side  of  this  equation  by  ()_  : 

ij 

* 

Qij  = / eiejda2n  ; (5_6) 

°2n 

we  must  prove  that  Qi;.  is  zero. 

In  fact,  if  follows  from  the  definition  (5-6)  that 
is  invariant  with  respect  to  an  interchange  of  the  two  axes 

and  x^  ; it  is  thus  the  same  regardless  of  whether  the  co- 
ordinate system  is  right-handed  or  left-handed.  Because  of  the 
spherical  symmetry,  is  invariant  with  respect  to  rotation 

and  to  reflection;  it  only  depends  on  the  geometrical  configura- 
tion. This  geometrical  conf i gurati on--2n+l  mutually  orthogonal 
axes--remai ns  unchanged  if  we  replace  the  x^-axis  by  its 
opposite  direction,  giving 


which  must,  therefore,  be  equal  to  Qi;j  . From  = - Qi;. 

we  get  = 0 and  hence  (5-5). 

The  reader  is  invited  to  make  this  reasoning  clear  to 
himself  for  the  case  of  three-dimensional  space  with  i = 1 
and  j = 2 . (In  this  case,  (5-5)  is  equivalent  to  the  ortho- 
gonality of  the  first-degree  harmonics.  Why?) 
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So  far,  we  have  restricted  our  considerations  to  the 
2n  + l anm  corresponding  to  the  same  degree  n . Let  us  now 
consider  two  different  degrees,  say  n and  n*  . Any  two  coef- 
ficients anm  belonging  to  two  different  degrees  will  be 
assumed  to  be  stochastically  independent.  Thus  the  probability 
space  a is  the  Cartesian  product  of  infinitely  many  unit 
spheres : 


O = o2  x o4  x „e  x „8  x ...  x „2n  x «2nt2  x ...  (5-7) 

where  o„  denotes  the  2n-dimensional  unit  sphere  in  R„  , . 
Thus  the  dimensionality  of  the  spheres  increases  with  in- 
creasing n , in  contrast  to  the  spherical  case  where  the 
probability  space  is  the  Cartesian  product  of  infinitely  many 
identical  circles. 

To  repeat:  in  the  present  model  any  two  different  anm 
are  uncorrel ated , but  for  different  reasons:  if  the  two  coef- 
ficients belong  to  different  degrees,  then  they  are  uncorrelated 
as  a consequence  of  the  statistical  independence;  if  the  two 
coefficients  belong  to  the  same  degree  n , then  they  are  un- 
correlated because  of  the  orthogonality  relation  (5-5). 

A second  model  of  an  ergodic  stochastic  process  is  ob- 
tained by  taking  the  probability  space  as  rotation  group 
space;  this  is  the  three-dimensional  analogue  of  the  second 
example  of  an  ergodic  process  considered  in  sec.  3.  In  view  of 
its  basic  importance  we  shall  devote  the  next  section  to  this 
model . 

6 . Rotation  Group  Space 

As  we  have  seen  in  sec.  3,  rotations  of  the  circle, 
which  constitute  the  rotation  group  in  two  dimensions,  are 
described  by  one  parameter  u ranging  from  0 to  2ir  : the 
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group  of  rotations  of  the  plane  forms  a one-dimensional  space, 
which  may  be  identified  with  the  unit  circle. 

Rotations  of  the  sphere,  which  make  up  the  rotation 
group  in  three  dimensions,  are  described  by  three  parameters, 
for  which  we  may  take  three  Eulerian  angles:  the  group  of 
rotations  of  three-dimensional  space  forms  itself  a three- 
dimensional  space,  whose  coordinates  are  the  three  Eulerian 
angles.  This  three-dimensional  rotation  group  space  cannot  be 
identified  with  the  unit  sphere.  This  is  in  contrast  to  the 
case  of  rotations  of  the  circle  and  accounts  for  the  greater 
complexity  of  the  present  case. 

Various  authors  use  various  definitions  of  Eulerian 
angles.  We  follow  the  definition  of  Synge  (1960,  p.18),  which 
is  fairly  widely  used  and  is  best  suited  for  the  present  purpose 
because  of  its  relation  to  the  spherical  coordinates  e,  x . 

Let  a rectangular  coordinate  system  XYZ  be  rotated 
into  a position  xyz  by  a general  spatial  rotation.  This 
rotation  is  split  up  into  three  successive  rotations  around 
coordinate  axes.  The  first  rotation  is  about  the  Z-axis  through 
an  angle  a ; it  transforms  the  XYZ-system  into  X ^Y  Z . The 
second  rotation  is  about  the  Y^axis  through  an  angle  0 ; thus 
we  obtain  a system  x2YiZ2  * we  rotate  about  the  Z2-axis 

about  an  angle  y in  a positive,  or  -y  in  the  negati ve, sense , to 
to  obtain  the  desired  system  xyz  . 

The  three  angles  A,  0,  y are  the  Euler  angles.  They 
may  be  illustrated  in  the  following  way  (Fig. 6).  The  angles  $ 
and  A are  the  usual  polar  coordinates  of  the  new  z-axis.  Let 
Q be  the  point  in  which  the  z-axis  intersects  the  unit  sphere, 
and  denote  by  7 and  y the  parallels  through  Q to  the  x- 
and  y-axis,  respectively.  Then  y is  the  angle  which  7 forms 
with  the  meridian,  positive  when  counted  counterclockwise  (the 
figure  shows  a negative  y ). 





Figure  6.  The  Euler  angles  A,  0,  *F 


In  the  usual  terminology,  the  angle 
else  than  the  azimuth,  counted  clockwise,  0 
x-di rection . We  shall,  therefore,  put 


and  consider  0,  A,  A as  our  final  Eulerian  angles 
These  three  angles  define  a point  w in  ro 
space,  which  we  shall  denote  by  n (we  shall  later 
it  as  our  probability  space);  we  thus  put 
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The  respective  ranges  are 


0 < 0 < IT  , 

0 < A < 2tt  , (6-3) 

0 s A < 2tt  . 

A rotation  defined  by  the  three  Euler  angles  (6-2)  will  be 
denoted  by  R . The  value 


* (0,  0,  0) 


(6-4) 


corresponds  to  the  identity  transformation  I,  leaving  the  axes 
XYZ  unchanged;  symbolically, 

Ro  = I . (6-5) 


A unit  vector  t defined  by  spherical  coordinates  e , 
A has  the  components 


t 


sine  cosa 
sine  sinx 
cose 


(6-6) 


It  will  be  symbolically  abbreviated  as 


t = (e.A)  ; (6-7) 

this  notation  has  already  been  used  before;  cf. equation  (4-11). 

The  rotation  (6-2)  transforms  the  vector  t into 
another  unit  vector,  which  we  shall  denote  by 


(6-8) 
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it  is  convenient  to  consider  R as  a rotation  matrix,  so  that 

(O 

(6-8)  is  the  usual  product  of  a matrix  and  a vector. 

Clearly,  the  Euler  angles  e,  A of  the  rotation 
are  completely  different  from  and  independent  of  the  coordinates 
0,  x of  the  vector  t . There  is,  however,  an  interesting 
relation  between  these  two  sets  of  quantities.  Form  the  triple 


t = (e,  a,  a)  , 


with  an  arbitrary  value  a between  0 and  2ir  , and 


(6-9) 


t = (-  0,-  A,-  a)  . 


Then  it  is  easily  seen  that 


(6-10) 


R.xt=  • 


(6-11) 


which  is  the  unit  vector  of  the  Z-axis.  Thus,  the  operation  R_t 
rotates  an  arbitrary  unit  vector  t into  the  Z-axis.  This 
simple  fact  will  be  of  importance  later  on. 

After  these  introductory  geometrical  considerations  we 
are  in  a position  to  construct  our  stochastic  process.  We  take 
a basic  function 


f(t)  = f ( 0 , A ) 


(6-12) 


on  the  unit  sphere  and  define  our  stochastic  process  by 


f (t,o>)  = f(Rut)  . 


(6-13) 


This  is  in  analogy  to  the  two-dimensional  case,  equations  (3-12) 
and  (3-14);  we  shall  also  follow  the  respective  developments 
in  sec.  3 as  closely  as  possible. 


Again,  the  functions  f(t,w)  differ  from  each  other 
only  by  a rotation  of  the  sphere;  they  are  not  essentially 
different.  Our  model  is  suited  to  represent  the  case  in  which 
we  have  only  one  realization  f(t)  but  wish  to  formally  use  the 
mathematical  techniques  of  stochastic  processes.  This  is  the 
case  of  the  terrestrial  gravitational  field,  where 

f(t)  = T(e,  x)  (6-14) 

is  the  anomalous  potential  at  sea  level.  The  choice  (6-13)  is 
intimately  connected  with  homogeneity  and  isotropy,  i.e.,  with 
invariance  of  essential  features  with  respect  to  rotations  . 
More  about  this  will  be  said  in  the  next  section. 

According  to  (6-13),  probability  space  n is  rotation 
group  space,  a point  wen  being  defined  by  the  three  Eulerian 
angles  (6-2).  The  expectation 


E{  *}  = /// ( - )dn  (6-15) 

(2 

is  an  integral  over  rotation  group  space.  The  integration  is 
to  be  extended  over  the  range  (6-3);  the  problem  is  to  find  a 
suitable  volume  element  dQ  , defining  a probability  measure. 

The  product  of  two  rotations  R is  again  a rotation. 

The  vector 

RjR2t  (6-16) 

is  obtained  by  rotating  the  vector  first  by  the  matrix  R2  and 
then  by  the  matrix  Rt  . Assume  that  a spherical  triangle  P1P2P3 
is  brought  by  a rotation  Rj  into  the  position  P]P2P3  » the 
configuration  (angles  and  sides)  of  both  triangles  is  obviously 
identical.  Let  t , t2 , t3  be  the  position  vectors  of  Pj,  P2 , 
P3;  all  are,  of  course,  unit  vectors.  Similarly  tj,  t2,  tj 
are  defined  (Fig. 7).  Then 


Figure  7.  Rotation  of  a configuration 


Each  of  the  vectors  t , t2 , t3  ; t^ , t2 , t^  can  be 

obtained  by  rotating  a fixed  unit  vector  e (for  which  we  may 

take,  for  instance,  the  unit  vector  of  the  X-axis)  by  a certain 

matrix  R (co  x ) , R(u>2),  ....  R(wj)  ; we  write  R(ui)  instead  of 

R to  avoid  two-level  subscripts.  Then,  for  i = 1,  2,  3, 

“i 


Combining  (6-17)  and  (6-18)  we  have 


Here  i = 1,  2,  3 , but  clearly  the  configuration  rotated  by 
Rj  can  have  any  number  of  points. 


Thus,  multiplying,  from  the  left,  a set  of  rotation 

matrices  R ( co ± ) , i = 1,  2,  3 by  a fixed  matrix  R1 

preserves  the  configuration.  The  geometrical  configuration  is 
invariant  with  respect  to  left  multiplication.  Similarly  we  may 
show  the  invariance  of  geometry  with  respect  to  right  multi- 
pi  ication. 

Homogeneity  and  isotropy  imply  that  the  essential  pro- 
perties depend  only  on  the  geometrical  configuration.  Therefore, 
also  the  probability  measure  must  be  invariant  with  respect  to 
right  and  left  multiplication.  It  can  be  shown  that  for  a compact 
group  such  as  the  rotation  group,  there  is  essentially  (apart 
from  a constant  factor)  only  one  group  measure  that  is  both 
right  and  left  invariant  (Smirnow,  III  1,  § 89).  Such  an  in- 
variant volume  element  in  rotation  group  space  is 

dV  = sinOdGdAdA  ; (6-20) 

this  will  be  proved  later  in  this  section.  The  total  volume 
of  group  space  is,  by  (6-3), 

2ir  it  2ir 

V = / / / sinededAdA  = 8*  , 

A=0  0=0  A=0 

so  that 

do  = sinededAdA  (6-21) 

8 it 

is  the  desired  element  of  probability  measure. 

Now  we  are  ready  to  attack  the  computation  of  expectations 
and  covariances.  The  expectation  E{f(t,u>)}  becomes 


t 


1 
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E{  f ( t ,o) ) } = J//f(Rut)dn 


fffURJ_Tt)da 


(6-22) 


in  view  of  right  invariance  and  using  (6-11).  However,  Ruez 
transforms  the  unit  vector  ez  of  the  Z-axis  into  the  unit 
vector  of  the  z-axis,  which  has  the  spherical  coordinates  0 
and  A (Fig. 6).  Hence, 


f(Ruez)  = f(0,A)  , 


(6-23) 


and  (6-22)  becomes 


. 2tt  it  2 it 

E{f(t,oi)>  = III  f (0,A)sined0dAdA  . 

8tt  A=0  0=0  A = 0 


We  integrate  over  A and  replace  ©,  A by  e,  A,  respectively; 
obviously,  the  symbols  for  the  integration  variables  is  irrelevant, 
The  result  is 


- 2n  it 

E{  f ( t ,w ) } = — l l f ( 6 , A ) s i n 6 d 6 d A , 


(6-24) 


A=0  0=0 


which  is  simply  the  average  over  f(t)  over  the  unit  sphere. 

It  is  zero  if  f(e,A)  contains  no  zero-degree  spherical  harmonic, 
which  corresponds  to  our  usual  assumption  that  the  stochastic 
process  under  consideration  is  centered. 
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Then,  by  (4-16),  the  covariance  function  becomes 


C(t,u)  = E{ f ( t ) f (u ) } 


(6-25) 


with 


t = (e,A)  , 


u = (e'.V) 


(6-26) 


More  explicitly  this  is  written 


C(t,u)  = //Jf(R  t)f(R  u)da  , 

n 


which,  because  of  right  invariance,  is  equal  to 


C(t,u)  = ///f(RMR.Tt)f(RuR_Tu)dn  . 


Now,  by  (6-11)  and  (6-23), 


(6-27) 


(6-28) 


f (Ru)R_Xt)  = f(0*A) 


(6-29) 


gives  the  value  of  f at  a point  P with  spherical  coordinates 
(0,A)  , whereas 


f(RwR_Tu)  = f (0  * ,A ' ) 


(6-30) 


denotes  the  value  of  f at  some  point  Q = (o', a')  situated 
at  the  spherical  distance  ^ from  P (Fig. 8).  That  the 
spherical  distance  ij>  between  P and  Q is  equal  to  the 
spherical  distance  between  the  points  t and  u as  given  by 
(4-23)  follows  from  the  invariance  of  the  configuration  with 
respect  to  the  rotation  R R . It  is  also  easily  seen  that 

O)  - T 


#1 


Figure  8.  Invariance  of  the  spherical  distance  \p 


if  a is  chosen  as  the  azimuth  from  t to  u 
(6-2)  will  be  the  azimuth  from  P to  Q . 

Thus  (6-28)  becomes 


A 1 ) si nedodAdA 


On  replacing  o,  A,  A by  e,  x 


x 1 ) si nededxda 
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Formally,  this  is  only  a change  in  the  symbols  for  the  inte- 
gration variables;  but  geometrically  the  meaning  is  now  altered 
profoundly:  a comparison  with  (4-25)  shows  that  this  is  simply 
the  average  M , so  that,  by  (4-26), 

C(t,u)  = Etf(t,«)f(u,io)}  = M{  f (t)f  (u) } = r(*)  ; (6-33) 


the  true  and  the  empirical  covariance  functions  are  identical. 

Exactly  as  in  the  two-dimensional  case,  the  reason  for 
this  identity  is  simply  that  probability  space  is  made  to  coin- 
cide with  rotation  group  space,  so  that  the  expectation  E and 
the  rotation  group  average  M are  identical;  our  present  model 
is  trivially  ergo die.  The  present  lengthy  discussion  has  been 
made  because  the  first  brief  presentation  of  this  model  in 
(Moritz,  1972,  1973,  sec.  9)  has  proved  to  be  too  sketchy. 

It  is  very  interesting  to  study  the  stochastic  behaviour 
of  the  spherical -harmonic  coefficients  a^nm(a))  in  the  present 
model;  the  situation  is  less  simple  than  in  the  two-dimensional 
analogue . 

Specializing  (4-10)  for  the  present  model  we  get  for 
fully  normalized  harmonics 


anm(“)  = 77  //f(t,u>)T?nm(t)do 
a 

■ -Jr  nfOjKJtx* 

a (6-34) 

= 1 // f(R  R t)TT  (R  t ) da 

4 it  J J ' id  -to  nm ' -u)  ' 
o 

= 77  Uf  (t)1?nm(R-ut)da  » 
o 


IL 


again  because  of  the  rotational  invariance  of  the  integral. 

It  is  now  well  known  how  spherical  harmonics  transform 
under  rotation.  References  to  the  numerous  literature  are  given 
in  (Aardom,  1969);  we  shall  mainly  follow  (Courant  and  Hilbert, 
1953,  pp, 535-545).  The  transformed  Legendre  harmonics  are 

simply  linear  combinations  of  all  2n+l  Legendre  harmonics  of 
the  same  degree  n : 

*»<'>-„*>  ’ l i <«-35> 

£.  = -n 

the  coefficients  Anm£  in  this  linear  combination  are  obviously 
functions  of  w , that  is,  of  the  rotation  parameters. 

On  substituting  (6-35)  and  integrating  we  get 

a (u)  = l A ,(«)a  0 , (6-36) 

nm'  ' nL  nmJi. v ' n£  ’ v ' 

l = -n 

where 

• 17  //m>*nt(t)d0  (6-37) 

o 

are  the  constant  coefficients  of  the  original  function  f(t)  . 

Eq.  (6-36)  is  the  three-dimensional  analogue  of  (3-24) 
and  (3-25).  Unfortunately,  the  explicit  expressions  of  ^njnA(“) 
are  rather  complicated.  So  we  shall  not  try  by  direct  computation 
to  verify  the  orthogonality  of  the  ^nm(“)  and  the  fact  that 
they  satisfy  (5-1):  this  follows  indirectly  from  the  identity 
between  the  empirical  and  the  true  covariance  function.  However, 
we  should  like  to  point  out  the  following  interesting  difference 
between  the  two-dimensional  and  the  spatial  case.  In  two 
dimensions,  the  stochastic  orthogonality  of  the  Fourier  coef- 
ficients has  been  a consequence  of  the  usual  orthogonality  of 


the  trigonometric  functions  over  the  unit  circle,  since  the 
transformation  coefficients  are  simply  given  by  cosku  and 
sinkw  . In  the  present  three-dimensional  model,  the  stochastic 
orthogonality  relations  (4-18)  are  not  a consequence  of  the 
orthogonality  relations  of  usual  spherical  harmonics  over  the 
unit  sphere,  as  might  be  expected  by  analogy.  In  fact,  rotation 
group  space  is  not  the  usual  unit  sphere,  but  is  described  by 
three  parameters. 

Now  it  is  very  curious,  however,  that  the  stochastic 
orthogonality  relations  in  three  dimensions  are,  in  fact,  a 
consequence  of  orthogonality  relations  for  spherical  harmonics, 
but  in  four-dimensional  space!  In  fact,  rotation  group  space 
may  be  identified  with  the  three-dimensional  "surface"  of  a 
unit  sphere  in  four-dimensional  space,  and  the  rotation  coef- 
ficients Anmjl(u>)  are  essentially  spherical  harmonics  in  this 
space;  see  below. 

We  shall  now  give  some  supplementary  information  on  the 
mathematical  structure  of  rotation  group  space.  The  reader  not 
interested  in  these  mathematical  details  may  proceed  directly 
to  the  next  section. 

Mathematical  complements.  According  to  (Smirnow, 
Vol.III/1,  sec. 90, p. 271) , the  invariant  integral  in  rotation 
group  space  has  the  form 


J//f(alt  a2,  a3)  j—  ^ da1da2da3 


/1-ara2’a3 

= III f(alt  a2,  a3)^-  dajda2da3  , 


(6-38) 


where  aQ , ax,  a2 , a3  denote  four  parameters  related  by 


2 2 2,2  , 

a„  + a.  + a0  + a,  = 1 

0 12  3 


(6-39) 


mm 
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These  four  parameters,  three  of  which  are  independent,  may  be 
used  to  describe  a rotation.  They  are  called  Eulerian  parameters 
(not  to  be  confused  with  Eulerian  angles)  and  are  related  to 
the  representation  of  rotations  by  quaternions. 

Geometrically  the  parameters  alt  a2,  a3,  a4  may  be 
interpreted  as  Cartesian  coordinates  in  a four-dimensional 
auxiliary  space;  then  (6-39)  describes  the  unit  sphere  in  this 
space . 

The  volume  element 

dV  = ~ da1da2da3  (6-40) 

O 

in  the  invariant  integral  (6-38)  is  nothing  else  than  the  three- 
dimensional  "surface"  element  of  this  four-dimensional  unit 
sphere.  In  fact,  the  ordinary  surface  element  of  the  unit  sphere 
in  three-dimensional  space,  expressed  in  terms  of  Cartesian 
coordinates  xyz  , is 

da  - * ^ dxdy  = ^ dxdy  (6-41) 

l/l-x^-y2 


(cf.  Smirnow,  V o 1 .II,  sec. 62,  p.176),  of  which  (6-40)  is  the 
four-dimensional  analogue. 

The  expression  of  the  Eulerian  parameters  a±  in  terms 
of  the  Eulerian  angles  A,  0,  v = -A  is  (Synge,  1960,  p.19, 
eq . ( 1 1 . 7 ) ) : 


a i = - sin|  sin^-^  , 


, _ „ . 0 . . A+A 

a 2 - si n-jj  co s- 


a 3 = c o s y si n^g— 


a = cos-*  cos-A_A 

O C 


(6-42) 


55 


I 


If  we  wish  to  express  the  group  volume  element  (6-40)  in  terms 
of  A,  0,  A we  have 

dV  = — J dedAdA  , 

ao 

where 

oa1  aa2  3a3 
J7T  ITT  JIT 

j - 

0 30  30  30 

3a1  3a2  3a3 
3 A 3 A 3 A 


is  the  Jacobian  determinant  of  the  transformation  (6-42). 
On  differentiation  of  (6-42)  we  find 


(6-43) 


(6-44) 


Me  subtract  the  first  row  in  (6-45)  from  the 
the  determinant  with  respect  to  the  third  row 
The  result  is 


a sin© 


Thus  (6-43)  becomes 


s i nededAdA 


P (cose)e 


It  is  slightly  different  from  (6-35)  but  is  essentially  equi 
valent  in  view  of 


cosrx  + i sinrx 


For  the  transformation  coefficients  S2'n  there  holds  ( ibid . , 
p.544) 


-A,r  _ -2nHn+«,,n+r/  . 

S2n  " V H2n  ^o*  ai*  a2*  a3'  ' 


(6-51) 


where 


2 2 . 2 . 2 
ao  + ai  + a2  + a3  ; 


(6-52) 


the  functions  H are  harmonic  polynomials  in  ak  (ibid. , p.542): 


32  4.  32  x 32  x 32  xHn+4,n+r  _ n 

~~2  + ~2  + 7"2  + ~ 2 ' "2n  " 0 

aa  3a.  3a0  3a, 

O 1 Z J 


(6-53) 


The  functions  S?;,r  are  defined  on  the  surface  of  the 
2n 

unit  sphere  (6-39)  and  form  an  orthogonal  system.  Courant  expresses 
them  in  terms  of  three  parameters  p,  a,  t,  which  are  related  to 
our  parameters  0,  A,  A by 


0 a-A 

T * 7 * p = ~T~ 


(6-54) 


7.  Statistical  Distributions  in  Rotation  Group  Space 


In  the  last  section  we  have  studied  rotation  group  space 
as  a probability  space  primarily  with  respect  to  covariances.  The 
covariance  theory  of  stochastic  process  is  what  is  needed  for 
linear  least  squares  prediction  and  estimation  problems;  it  can 
be  treated  without  explicit  reference  to  the  underlying  statistical 
distributions,  of  which  only  the  moments  of  first  order  (mean 
values),  and  of  second  order  (variances  and  covariances)  are 
needed . 
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Even  in  least-squares  prediction  and  collocation,  how- 
ever, the  distributions  of  relevant  quantities  are  required  if 
we  wish  to  perform  statistical  tests.  A1 ready  to  answer  very 
elementary  but  meaningful  questions  we  need  distributions. 

Such  a question  is,  for  instance:  What  is  the  average 
global  frequency  of  a l°x  1°  mean  gravity  anomaly  situated 
between  - 28  and  - 36  mgal?  This  question  is  answered  by  the 
histogram  of  Fig.  9;  the  frequency  is  the  number  of  l°x  1° 
anomalies  having  magnitude  within  a specified  interval,  divided 
by  the  total  number  of  l°x  1°  anomalies.  Clearly,  such  a fre- 
quency can  be  considered  as  a measure  of  the  probability  that 
a l°x  1°  mean  anomaly  lies  between  - 28  and  - 36  mgal. 

A similarly  meaningful  question  would  be:  What  is  the 
probability  that  a l°x  1°  mean  Ag-value  lies  between  - 28  and 
- 36  mgal  and  that  the  mean  value  of  the  geoidal  height  N for 
the  same  l°x  1°  block  lies  between  25  and  30  meters? 


Figure  9.  Number  of  l°x  1°  mean  gravity  anomalies  having  magnitude 
within  specified  interval.  After  (Rapp , 1977  ,p . 5) . 
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To  answer  such  and  related  question,  we  must  construct 
appropriate  distribution  functions  for  Ag  , for  Ag  and  N 
jointly,  etc.  The  distribution  density  for  Ag  will  be  a 
continuous  analogue  of  the  histogram  of  Fig. 9.  To  find  an 
appropriate  probability  space,  let  us  note  that  the  number  of 
relevant  l°x  1°  mean  anomalies  is  counted  regardless  of  the 
position.,  jn  the  earth's  surface,  of  the  l°x  1°  blocks  under 
consideration.  All  positions  on  the  sphere  are  treated  equally: 
again  we  have  homogeneity  and  isotropy.  Thus,  rotation  group 
space  is  seen  to  be  the  proper  probability  space  also  as  a basis 
for  the  mathematical  description  of  statistical  distributions. 

The  nature  of  a statistical  distribution  is  best  illus- 
trated by  the  case  of  a function  of  one  variable.  Therefore,  we 
shall  again  start  with  the  group  of  rotations  of  the  circle  (the 
two-dimensional  rotation  group),  which  is  parametrized  by  one 
variable  to,  0 < to  < 2-rr  . Probability  space  n is,  therefore, 
the  unit  circle,  and  the  element  of  probability  measure  is  ^ dw, 
which  is  clearly  left  and  right  invariant.  Let  the  random 
function  under  consideration  be  denoted  by  f(oj)  . 

We  plot  a)  along  the  horizontal  axis  of  a graph;  then 
f(w)  is  defined  for  0 ^ <o  < 2 it  (it  could,  of  course,  be  con- 
tinued periodically  for  other  abscissas). 

Then  the  distribution  function  $(x)  is  defined  by 

4>(x)  = Prob{f(a>)  < x)  (7-1) 

as  the  probability  that  f(a>)  takes  a value  smaller  than  x . 

It  is  the  measure  of  all  values  of  to  for  which  f(u>)  < x ; 
this  measure  is  obviously  a function  of  x . In  the  situation 
shown  in  Fig.  10,  f ( c*> ) < x if  to  is  contained  in  the  interval 
AB  or  in  the  interval  CD  ; thus 

Prob(f(<o)  < xl  = + TTD)  , 


i 


i 
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Figure  10.  Distribution  of  a random  function  f(w) 


denoting  the  length  of  AB  and  the  factor  1/2*  serving 
to  make  the  measure  of  the  total  interval  from  0 to  2*  equal 
to  unity. 

Generally  we  may  write 

*(*)  = 77  / d“  • (7-2) 

f ( U)  ) < X 

the  integral  being  extended  over  those  points  u>  for  which 
f(u>)  < x . 

The  derivative  of  (7-1)  with  respect  to  x gives  the 
probabil ity  density 


-3 


which  has  an  even  more  intuitive  geometrical  meaning.  For  a 
differential  dx  (operations  with  differentials  are  justified 
in  the  usual  way)  we  have 


d> ( x ) dx  = d$(x)  = Probix  < f(w)  < x + dx)  . 


(7-4) 


According  to  Fig.  11,  which  represents  the  same  function,  f(u>) 
assumes  a value  between  x and  x + dx  if  x lies  in  one 
of  the  small  intervals  A'A,  BB',  C'C,  or  DD1,  so  that 


d$(x)  = i-(  A' A + BE'  + C'C  + T5TJ' ) , 


(7-5) 


*(*)  = + bf*  + rr  + M' ) , 


(7-6) 


D D’  2* 


Figure  11.  Definition  of  the  distribution  density 


which  gives  an  intuitive  geometrical  interpretation  of  the 
distribution  density  4>(x)  . 

The  distribution  function  $(x)  is  a monotone  non- 
negative function  of  x,-“><x<«,as  shown  in  Fig.  12.  It 
is  identically  zero  for  x < frain  (there  is  no  w for  which 
f (to ) < f . ) and  identically  one  for  x > f (for  al  1 u. 

there  is  f(u>)  < x if  x > f ) . 

The  statistical  expectation  E of  the  random  variable 
f can  now  be  computed  in  two  ways:  by  means  of  the  probability 
measure  du>/27r  : 


E { f } = Jf (w)dm  , (7-7) 

or  by  means  of  the  distribution  function 

oo 

E ( f > = /xd$(x)  . 


(7-8) 
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Geometrically,  both  expressions  give  the  area  under  the  curve 
f(«)  in  Fig. 10;  (7-7)  corresponds  to  the  Riemann  partitioning 
and  (7-8)  corresponds  to  the  Lebesgue  partitioning  of  the  same 
definite  integral;  therefore  (7-7)  and  (7-8)  are  identical; 
cf.  (Feller,  1966,  p. 115-116)  and  (Kolmogorov  and  Fomin,  1970, 
p . 293 ) . 

Let  us,  finally,  consider  the  stochastic  process  (3-12), 

f(t.«)  = f(t+u>)  . (7-9) 

In  view  of  the  rotational  invariance,  the  distribution  function 
of  f ( t+oi ) , for  fixed  t , is  the  same  for  all  t , hence,  is 
the  same  as  for  t = 0 : 

Prob{f(t+u>)  < x}  = Prob(f(u)  < x}  . (7-10) 

This  is  immediately  seen  from  Fig.  11:  replacing  t by  t+w 
means  only  a translation  of  the  figure  as  a whole  to  the  right 
or  left;  the  length  of  the  intervals  A 1 A , ...  and  hence  (7-5) 
or  (7-6),  remaining  unchanged. 

What  is  more,  we  may  also  write 

$(x)  = Measif (t)  < x}  , (7-11) 

using  only  the  sample  function  f(t)  defined  on  the  "space 
circle"  0 < t < 2n  with  measure  "Meas"  defined  by  its  element 
dt/2n  , without  any  probabilistic  i nterpretati on ; this  follows 
immediately  by  replacing  w in  (7-1)  by  t . This  is  formally 
very  simple,  but  conceptually  of  fundamental  importance:  it 
shows  that  we  may  consistently  work  with  one  basic  sample  func- 
tion f(t)  only  and  still  avail  ourselves  of  the  formal  advan- 
tages of  probability  theory. 
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Distributions  in  three-dimensional  rotation  group  space, 


After  these  preliminaries  we  come  to  the  geodeti cal ly  relevant 
case  of  three-dimensional  rotations.  The  basic  ideas  remain  the 
same,  though  the  notation  is  more  cumbersome. 

Let  f(a>)  be  a real-valued  random  function;  the 
argument  is  defined  by  (6-2).  Then  the  distribution  function 
t(x)  of  f is 


4>(x)  = Prob{f(a>)  < x}  . 


(7-12) 


It  should  be  noted  that  x is  a one-dimensional  real  variable, 
- oo  < x < 00  , though  ui  denotes  a point  in  three-dimensional 
rotation  group  space  n . 

Consider  now  the  random  function  (6-13), 


f(t,u>)  = f(Rwt)  with  t = ( e , X ) . 


(7-13) 


In  view  of  the  rotational  invariance,  the  distribution  function 


*(x)  = Prob{f(Rajt)  < x} 


(7-14) 


does  not  depend  on  t . Following  the  reasoning  that  leads  from 
(6-22)  to  (6-24)  we  find  that 


<t»(x)  = Meas(f(e,x)  < x> 


(7-15) 


The  measure  "Meas"  is  surface  measure  on  the  unit  sphere, 
normalized  by  the  factor  l/4ir  ; its  element  is,  as  usual. 


— do  = I—  sinededx  . 
4w  4 it 


(7-16) 


Just  as  in  (6-24),  there  is  no  longer  an  explicit  dependence  on 
the  azimuth  variable  A . 

For  the  distribution  density 


♦(x)  - *'(x) 


(7-17) 


we  have  again  a geometrical  interpretation  (Fig. 13).  Draw  the 
neighboring  contour  lines 


f (0 ,x)  = x = const.  , 
f(e,A)  = x + dx  = const. 

on  the  sphere;  they  will,  in  general,  consist  of  several  un- 
connected closed  curves.  Let  the  areas  between  these  neighboring 
closed  lines  be  denoted  by  dA2,  dA3,  ...  (hatched  in 

Fig. 13).  Then 


f-x+dx 


f-x+dx 


f»X+dX, 


Figure  13.  Geometrical  interpretation  of  the 
distribution  density 
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4>(x)dx  = + dA2  + dA3  + ...)  . 


(7-18) 


The  distribution  function  4>(x)  itself  can  be  expressed 
in  a form  analogous  to  (7-2): 


# ( x)  =4--  //  sinededA  . 

* f (0,  A)  <K 


(7-19) 


Another  basic  problem  is  the  determination  of  the  joint 
distribution  of  two  functions  f and  g on  the  sphere,  say,  of 


f(e,x)  = T(e,x)  , 
g(e,x)  = Ag(e,x)  , 


(7-20) 


T and  Ag  denoting  the  disturbing  potential  and  the  gravity 
anomaly,  respectively.  The  joint  distribution  function  is 


*(x,y)  = Measif (e,x)  < x,  g(e,x)  < y> 


(7-21) 


The  corresponding  density  is 


♦ (x.y)  = ; 

3x3y  * 


(7-22) 


it  may  be  geometrically  illustrated  as  follows.  Draw  the  contour 


as  well  as  the  contour  lines 


g(0.x)  ■ y , 
g(e,x)  = y + dy 


(7-24) 


(Fig.  14).  The  ribbons  formed  in  this  way  intersect  in  areas 
dAj,  dA2,  dA3 , ...  (hatched  in  Fig.  14),  and 


<f>  (x,y)dxdy  = ^(dAj  + dA2  + dA3  + ...)  . (7-25) 


A final  example  will  indicate  how  an  azimuth-dependent 
situation  can  be  handled.  Consider  the  problem  of  the  joint 
distribution  of  gravity  anomalies  at  two  points  that  are  at  a 
spherical  distance  * apart: 

*(x,y)  * Prob{f(t,w)  < x,  f(u,u>)  < y}  (7-26) 


Figure  14.  Joint  distributions 


~m r 


jtj  ^..i.-L--crr*r 
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where 


ill  = angle(t.u)  = const. 


(7-27) 


We  have 


t = (e,x)  , (7-28) 

u = (0  * , A * ) , (7-29) 

where  the  condition  (7-27)  can  be  written  in  the  form 

cosecose'  + si nesine ' cos (x ' -x ) = cosij»  = const.  (7-30) 

Then  (7-26)  can  be  expressed  as  the  integral 

<J>(x,y)  = — ~2  ///  sinededxda  , (7-31) 

8tt  B (x  , y ) 


where  the  integration  is  extended  over  the  region  B(x,y) 
defined  by  the  inequalities 

f(9»x)<x,  (7-3 

9(e,,x’)  < y ; 

0',  x'  are  expressed  as  functions  of  (e,X,a)  by  the  tri- 
gonometric relations 


cose'  = cosecose  + sinesin^cosa  , 
sin(x'-x)  = sinijisina/sine'  , 


(7-33) 


which  follow  from  the  spherical  triangle  of  Fig.  15.  The  integral 
(7-31)  is  analogous  to  (7-2)  and  (7-19);  the  probability  measure 
sinedodAdA  has  been  replaced  by  ’’spatial"  (surface  plus  azimuth) 


Figure  15.  The  basic  spherical  triangle 


measure  sinededxda  in  the  same  way  as  (6-31)  has  been  replaced 
by  (6-32). 

These  three  examples  illustrate  the  basic  principles  of 
the  determination  of  single  and  joint  distributions.  Other  cases 
can  be  handled  similarly.  In  any  case,  we  can  operate  with 
"spatial"  functions  f(8,x),  g(e,x),  ...  only. 

In  practice,  the  functions  f(e,x)  are  usually  re- 
presented by  discrete  mean  values  (say,  5'x  5'  or  l°x  1°  block 
averages),  and  the  integrations  are  to  be  replaced  by  sums. 


8.  The  Meaning  of  Statistics  in  Collocation 

Are  gravity  anomalies  a stochastic  phenomenon?  There  are 
different  answers  to  this  question. 

To  get  one  answer,  consider  gravity  at  one  observation 
station  and  observe  it  repeatedly.  Assume  that  the  measuring 

I ( 
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errors  are  negligibly  small,  and  remove  known  geophysical  effects, 
especially  tidal  ones.  Then  the  results  of  measurements  at 
different  times  will  be  practically  constant.  We  conclude:  gravity 
is  a deterministic,  not  a stochastic  phenomenon. 

This  way  of  looking  at  the  problem,  repeating  the  same 
experiment  under  identical  conditions,  is  the  way  we  look  at 
random  measuring  errors  and  at  many  other  "stochastic"  physical 
phenomena,  the  classical  case  being  the  repeated  throw  of  a die. 

If  the  experimental  results  vary  randomly,  then  we  have  a 
genuinely  stochastic  phenomenon.  Under  the  assumption  that  the 
outcomes  of  the  repeated  experiment  are  independent  of  one  another, 
we  have  the  scheme  of  repeated  trials,  fundamental  in  probability 
theory . 

There  is,  however,  also  another  way  of  looking  at  the 
question  of  stochasticity  of  gravity  anomalies.  They  are  caused 
by  mass  anomalies,  visible  and  invisible  ones.  These  mass 
anomalies  show  some  regular  features,  for  instance,  mountain 
chains  extending  in  a regular  fashion  from  north  to  south.  After 
removing  known  irregularities,  however,  the  residuals  are  rather 
irregular.  It  is  difficult  to  recognize  a regular  pattern.  We  may 
say,  with  some  justification,  that  the  residual  gravity  field 
is  caused  by  randomly  distributed  mass  anomalies.  In  this  sense, 
gravity  anomalies  (after  subtraction  of  known  regular  trends) 
may  be  said  to  be  random,  perhaps  even  stochasti cal  . The  ran- 
domness exi sts  here  not  with  respect  to  time,  as  it  was  in  the 
first  case  (measurements  at  the  same  point  but  at  different 
times),  but  with  respect  to  space  (measurements  at  the  same 
time  but  at  different  points).  The  random  behavior  is  more  or  less 
independent  of  position  on  the  sphere  and  of  direction:  it  is 
homogeneous  and  isotropic. 

Thus,  the  global  anomalous  gravitational  field  may  be 
irregular  enough  to  be  considered  as  a realization  (a  sample 
function)  of  a stochastic  process.  Is  this  sufficient  for  saying 
that  " the  anomalous  gravitational  field  is  a stochastic  process"? 


Personally,  I do  not  think  so,  because  there  is  no  other  physical 
realization:  there  is  only  one  Earth. 

The  situation  might  well  be  compared  the  problem  of  a 
global  statistics  of  the  human  population.  There  is  a temporal 
variation,  but  it  is  systematic  (expansion)  rather  than  random, 
and,  not  being  associated  with  the  Club  of  Rome,  I shall  not 
consider  it  here.  But  we  have  random  variations  from  one  human 
individual  to  the  other.  There  are  regular  trends--col or  of  skin, 
political  and  religious  beliefs--but  there  are  genuinely  irregular 
features  left,  distributed  over  the  human  population  and  hence 
over  the  earth's  surface.  This  is  not  completely  unlike  the  sur- 
face distribution  of  gravity  anomalies  (although  the  analogy, 
if  pushed  too  hard,  quickly  becomes  nonsense). 

Is  it  permitted  to  study  the  global  population  statistics 
at  a given  time,  to  calculate  various  statistical  distributions? 
Every  one  will  answer  this  question  in  the  affirmative,  although 
there  is  only  one  global  population.  All  statistical  distri- 
butions are  simply  calculated  on  the  basis  of  this  population. 

I think  the  statistics  of  the  gravitational  field  must 
be  handled  similarly.  We  simply  must  take  seriously  the  fact 
that  there  is  only  one  field,  and  compute  the  whole  statistics 
from  this  one  field  only. 

The  appropriate  mathematical  apparatus  for  studying  the 
"second-order  statistics"  (variances  and  covariances)  of  the 
gravitational  field  is  thus  Norbert  Wiener's  (1930)  "covariance 
analysis  of  individual  functions"  (Doob,  1949,  sec.l).  This 
model  is  implicit  in  almost  all  geodetic  work  in  this  field  (Kaula, 
1959,  1967  ; Heiskanen  and  Moritz,  1967)*,  explicitly  it  was 
formulated  in  (Moritz,  1973,  sec. 8).  It  essentially  uses  the  idea 
of  homogeneity  and  isotropy.  For  the  sphere,  homogeneity  and 
isotropy  really  form  a single  compound  notion,  namely,  in- 
variance under  rotations  (in  contrast  to  the  plane,  where 
homogeneity,  invariance  under  translation,  and  isotropy,  in- 
variance under  rotation,  are  separate  notions!);  this  motivates 
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the  introduction  of  the  three-dimensional  rotation  group. 

The  present  report  attempts  to  extend  Wiener's  idea  * 

beyond  a second-order  theory,  in  such  a way  as  to  obtain  a 

I 

complete  statistical  theory  including  statistical  distributions. 

This  has  been  done  in  the  last  section;  it  has  been  seen  that 
also  distributions  can,  in  fact,  be  obtained  from  one  given 
function  only. 

Formally,  this  theory  can  also  be  interpreted  within 
the  framework  of  stochastic  processes,  as  our  ergodic  Second 
Model.  This  is,  of  course,  independent  of  the  question  whether 
the  anomalous  gravitational  field  is  "really"  a stochastic 
process  in  some  physical  sense.  Probability  theory  simply  serves 
to  provide  a convenient  mathematical  formalism.  In  this  sense, 
formal  probabilistic  techniques  have  been  successfully  applied 
not  only  in  analytical  mechanics  with  a large  number  of  particles 
(Khinchin,  1949),  but  even  in  analysis  and  number  theory  (Kac, 

1959a, b) . 

The  problem  whether  and  in  which  respect  the  anomalous 
gravitational  field  is  a "genuinely  stochastic  phenomenon"  will 
be  answered  di fferently  by  different  people,  depending  on  their 
scientific  outlooks.  Even  with  respect  to  the  philosophical 
meaning  of  "probability"  and  "stochastic  phenomenon"  there  are 
many  different,  even  quite  opposite,  opinions,  as  is  seen  by 
comparing  books  such  as  (Gnedenko,  1967,  and  (de  Finetti,  1972); 
still,  the  mathematical  formalism  is  the  same. 

Similarly,  the  mathematical  formalism,  proposed  here 
as  a statistical  background  of  collocation,  is  independent  of 
how  serious  we  take  the  stochastic  character  of  the  gravitational 
field.  Even  if  we  rigorously  deny  this  stochasticity,  we  can 
still  accept  the  formal  statistical  analysis  presented  here:  we 
then  have  "statistics  without  stochastics". 

The  statistics  of  measuring  errors  can  be  incorporated 
without  problems,  in  the  way  described  in  (Moritz,  1973,  secs. 8 
and  9):  the  combined  phase  space  is  the  Cartesian  product  of 
rotation  group  space  and  of  the  probability  space  of  the 
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measuring  errors;  a combined  distribution  function  is  simply 
the  product  of  the  distribution  function  of  the  field  quantities 
under  consideration  and  the  distribution  function  of  their 
measuring  errors;  and  the  averaging  operator  to  be  used  is  the 
(commutative)  product  of  the  rotation  group  average  M and  of  the 
statistical  expection  E referring  to  the  probability  space  of 
the  measuring  errors.  If  we  limit  ourselves  to  a second-order 
theory,  then  the  approach  of  Sansd  (1978)  is  logically  particu- 
larly satisfactory. 

If  the  approach  of  sections  6 and  7 is  accepted,  then  a 
detailed  theory  of  statistical  distributions  for  geodetically 
relevant  quantities,  such  as  gravity  anomalies,  geoidal  heights, 
and  deflections  of  the  vertical,  could  be  developed  and  applied 
to  the  statistical  testing  of  the  results  of  least-squares 
col  1 ocation . 
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